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I .  INTRODUCTION 


Considered  here  is  the  theoretical  modeling  of  turbulent 
mixing  at  density  discontinuities  in  nonsteady  compressible 
flows.  It  is  well  known  that  if  a  density  discontinuity  (or 
a  strong  density  gradient)  occurs  in  a  pressure  gradient  of 
the  opposite  sign,  then  the  flow  field  is  hydrodynamically 
unstable  in  the  Rayleigh-Taylor  sense  (Taylor,  1950).!  Small 
perturbations  occurring  in  such  a  region  will  amplify  with 
time,  and  if  the  perturbations  become  large  enough,  they  can 
lead  to  a  local  breakdown  in  the  well-ordered  flow;  i.e., 
they  can  lead  to  turbulence. 

One  of  the  more  interesting  examples  of  density  gradients 
working  against  pressure  gradients  to  cause  instabilities 
occurs  in  blast  waves  driven  by  solid  high  explosives  (HE).* 
After  the  detonation  wave  breaks  out  of  the  charge,  the  HE 
combustion  gases  expand  to  a  high  velocity  (^6  km/s) ,  pushing 
an  air  shock  ahead.  From  one-dimensional  inviscid  calcula¬ 
tions  of  this  problem  (Brode,  1957), 2  we  know  that  a  positive 
pressure  gradient  is  formed  throughout  the  flow  field.  Such 
calculations  also  indicate  that  there  is  a  large  density  jump 
{pH£/pair  can  be  as  large  as  70)  across  the  contact  surface. 
However,  from  high-speed  photography  of  HE  experiments,  we 
know  that  this  contact  surface  (which  theoretically  should  be 
smooth)  actually  develops  an  irregular  shape  indicating  the 
growth  of  instabilities.  Evidence  of  mixing  at  the  contact 
surface  can  be  inferred  from  test  results  which  show  that  the 
HE  gases  react  with  the  shocked  ambient  gases  and  release 


* 

Similar  gradients  can  occur  in  "shock- tube- type"  flows  driven 
by  gases  at  very  large  initial  densities. 

1.  Sir  G.  I.  Taylor  (1950),  "Instability  of  Liquid  Surfaces 
when  Accelerated  in  a  Direction  Perpendicular  to  Their 
Plane,"  Proceedings  of  the  Royal  Society  (London),  A*^l, 
pp.  192-196. 

2.  H.  L.  Brode  (1957),  A  Calculation  of  the  Blast  Wave  from 
a  Spherical  Charge  of  TNT,  Rand  Corporation,  RM-1965. 


heat  (i.e.,  afterburn)  if  the  ambient  gases  contain  oxygen 
(Matle,  1959,  and  Filler,  1956)^»4.  Since  this  heat  release 
can  be  of  the  same  order  as  the  energy  released  by  the  detona¬ 
tion  wave  (e.g.,  for  TNT  postdetonation  energy  release  is 
about  2.5  times  the  detonation  energy),  it  will  affect  the 
blast  wave  flow  field  and  must  be  taken  into  account  for 
accurate  numerical  simulations  of  such  flows. 

This  report  describes  the  numerical  work  that  RDA  has 
performed  to  simulate  the  kinematics*  of  the  turbulent  mixing 
associated  with  such  HE-driven  blast  waves.  A  k-e  turbulence 
model  was  programmed  and  combined  with  RDA's  minimal-diffusion 
Flux-Corrected  Transport  (FCT)  module  which  solves  the  gas 
dynamic  conservation  laws  on  a  sliding  grid.  Source  terms 
were  included  in  the  k-e  equations  to  model  the  above- 
described  Rayleigh-Taylor  mechanism  to  generate  turbulence. 

The  governing  equations  are  described  in  Section  II.  Prelim¬ 
inary  numerical  studies  (as  described  in  Section  III)  indicate 
that  this  model  does,  indeed,  generate  turbulent  kinetic 
energy  at  the  HF./air  interface  of  HE-driven  blast  waves. 
Conclusions  and  recommendations  are  offered  in  Section  IV. 


Energetics  (i.e.,  heat  release)  was  neglected  for  the 
preliminary  work. 

3.  C.  C.  Matle  (1959),  The  Contribution  of  Afterburning  to 
the  Air  Blast  from  Explosives,  NAVORD  Report  6234. 

4.  W.  S.  Filler  (.1956),  "Post-Detonation  Pressure  and  Thermal 
Studies  of  Solid  Explosives  in  a  Closed  Chamber,"  Sixth 
(Int)  Symposium  on  Combustion,  pp.  648-657. 
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II.  FORMULATION 


The  time  evolution  of  the  mean  flow  is  governed  by  the 
conservation  of  mass,  momentum  and  energy  which  may  be 
written  as  follows  for  point  ( j=2) ,  line  ( j=>l ) ,  or  plane 
(.1=0)  symmetric  flow: 


—  p  +  -if  (r^pu)  =  0 
3t  w  i  3r  v  H 
r 

W  “>u)  +  Tr  *  -Jr  +  I  pk) 

K  w  &  (r3‘"lE)  ■  5?  t3“(p  +  !  o'1) 


i„  /_J_  It 

T\°t.h  3r 


1  3k\ 

°t,k  3r/ 


+  Q 


(1) 

(2) 

(3) 


where  mass-averaged  properties  of  the  mixture  are 
p  o  density 

u  c  velocity  ,  2 

E  “  total  energy  *  e  +  ^  u  +  k 
e  «  internal  energy 
p  ®  pressure 

U  *  rate  of  energy  release  per  unit  volume 
h  *  enthalpy 

and  the  turbulenco  parameters  are 

k  a  turbulence  kinetic  energy  per  unit  mass 
°T,fi  *  turbulent  Prandtl  number  for  the  variable  8 
UY  *  total  viscosity  =  +  C*,  pk2/e 

«  molecular  viscosity 

For  a  mixture  of  N  component  gases,  one  must  solve  N-l  species 
transport  equations  in  addition  to  the  mass  conservation  law 
for  ho  mixture  (Eq.  1).  The  species  transport  equations  have 
the  form 


3 

3r 


(4) 
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where  Za  -  mass-averaged,  mass  fraction  of  component  a 

*  p  /p  (while  I  f  =1) 
u  a  cx 

*  • 

SA  *»  rate  of  creation  of  component  a  (with  Z  S  »  0) 
cx  ci  a 

c 

As  shown  by  Issa  (1980),  a  one-dimensional  k-e  turbulence 
model  can  be  used  to  simulate  turbulent  mixing  driven  by  the 
Rayleigh-Taylor  mechanism.  The  transport  equations  for  the 
turbulence  kinetic  energy,  k,  and  the  turbulence  kinetic 
energy  dissipation  rate,  e,  acquire  the  following  form: 


pe  (5) 


Ti  !ek)  *  7J  jjU^oku)  =  +  G  - 

i  £  (r3peu)  -  i  ±  (r2  ~  ||)  f  CXG  /k  -  C2B£2/k 


where 


G  «  -1  — _ _  _  n  ^  5p  3 


‘nil?  (r3u)  **  C3 


TP  T?  ^ 


(6) 


(7) 


Here  G  represents  a  source/ sink  for  turbulence  kinetic  energy. 
The  first  term  in  G  is  proportional  to  the  flow  divergence 
and  tends  to  kill  turbulence  in  expanding  flows.  The  second 
term  is  proportional  to  the  product  of  the  moan  pressure  and 
density  gradients.  When  those  gradients  have  opposite  signs, 
then  G  can  become  positive.  As  will  be  shown,  the  pressure 
and  density  gradients  do  have  such  opposite  signs  near  the 
contact  surface  of  an  HR-driven  blast  wave  and  turbulence 
indeed  grows  there,  in  other  regions  of  the  flow  (e.g.,  at 
the  shocks  and  in  the  rarefaction  wave),  the  pressure  and 
density  gradients  have  similar  signs  and  turbulence  is 
suppressed. 


5.  R.  l.  tssa  (1980),  Modeling  of  Turbulent  Mixing  at 

Density  Discontinuities  i"n  Nonsteady  Compressible  Flows , 
R  *  D  Associates.  RDA-TR-107605-001. 
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III.  NUMERICAL  STUDIES 


The  k-e  turbulence  equations  described  in  the  previous 
section  were  programmed  on  the  computer  and  incorporated 
into  the  one-dimensional  hydrocode,  RDAFCT.  This  cods  uses 
the  FCT  module  labeled  JPBFCT^  to  solve  general  conservation 
laws  on  a  sliding  grid: 


3 

at 


6A(t) 


w(u-u  ) 
5 


5A(t) 


dA  +  <ji  \ 


dA 


6A(t) 


(8) 


where  w  «  (p,  j  j,  pE,  pf,  etc.)  and  the  grid  motion  satisfies 
dA  (t)  =  Ug(t)dA(t).  The  module  is  called  success'' vely  for  each 
conservation  equation  (i.o.,  Eos.  (1)  to  (6))  to  advance  the 
flow  field  one  time  step.  NRL's  newest  transport,  diffusion 
and  antidiffusion  coefficients  are  employed  to  give  sixth-order 
phase  and  fourth-order  amplitude  properties  for  the  numerical 
algorithm. 


a .  Initial  cond i tions — As  a  test  problem,  we  chose  to 
simulate  the  evolution  of  the  flow  field  associated  with  the 
detonation  of  a  l~lt>  spherical  charge  of  FRX-94Q4  in  air. 
Predetonation  conditions  for  the  charge  were 

Rc  »  1.88973  cm 

a  *  1.84  q/ctn3 
o 

while  the  ambient  air  conditions  were  taken  as 

P  «  1.01325  x  1 0°  dy/cro2 

-  1.22s?  x  10" 5  g/cm3 

T  =•=  288. 4°K 
o 

It  was  assumed  that  the  charge  was  center-detonated .  Initial 
conditions  for  the  numerical  calculations  were  taken  at  the 
instant  the  detonation  wave  reached  the  charge  radius.  It 
wad  assumed  that  the  flow  field  corresponded  to  that  of  an 
ideal.  Chapman- Jouquot  {CJ)  detonation  wave:  this  field  was 
obtained  from  subroutine  CMDET  which  provides  the  similarity 
solution  for  any  given  explosive  (Kuhl,  19?8).6 7 


6 .  JT~F7~Uor is  {1976),  F lux-Corrected  Transport  Modules  for 
Solving  General iged  Continuity  Equations,  Naval  Research 
Laboratory,  Memorandum  Report  3237’. 

7.  A.  L.  Kuhl,  M.  U.  Seizew  (197R),  Analysis  of  Ideal,  Strong, 
Chapman- Jouguet  Detonations ,  TRW  Report  78.4735.9-13. 
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b.  Equation-of-state — The  detonation  products  were 
modeled  by  the  Jones-Wilkins-Lee  (JT4L)  equation-of-state 
(Dobratz,  1974)®  which  provides  the  pressure  as  a  function 
of  density  and  internal  energy  in  the  form 

p(p,e)  =  a/i-  R^W<-Vo/P> 


+  B(1'V;)eXP('R2Po/p)  +  “pe 

The  JWL  parameters  for  PBX-9404  are 

A  =  8.545  Mbars 
B  =  0.20493  Mbars 
%  =  4.60 

R2  =  1.35 
to  =  y-1  =  0.25 


while  the  CJ  detonation  parameters  for  this  explosive  are 

PCJ  -  po  “cj2/(r+1)  ■  370  kbaf 
PCJ  =  PQ  (r+l)/r  «  2.485  g/cm-3 

W__  =  8.8  km/s 
CJ 


=  w__/(r+l)  =2.28  km/s 

CJ  -t  ^ 

=  E0/P0  =  5.543  x  10-  erg/g 


FCJ  =  2‘85 

ecj  =  qcj  +  1/2  ucj2  *  8-142  x  1C)10  erg/g 


Local  thermodynamic  equilibrium  was  assumed  for  the  air.  The 
pressure  was  related  to  the  density  and  internal  energy  by  an 
effective  gamma: 

p  =  (Y  -i)  Pe  (10) 
6 


where  Ye  =  Ye(P»e)  was  obtained  from  a  table-lookup  subroutine 
which  uses  Gilmore's  data  for  real  air  (Gilmore,  1955).9  Two 
gas  species  were  considered:  air  and  detonation  products. 
Computational  cells  containing  air  were  initialized  with  f=0, 
while  cells  containing  detonation  products  were  initialized 
with  f=l.  A  transport  equation  was  then  included  for  the 
species  fraction  f.  In  the  present  nonreactive  case,  f  served 


8T  B.  Dobratz  (1974),  Properties  of  Chemical  Explosives 
and  Explosive  Simulants,  Lawrence  Livermore  Laboratory 
Report  UCRL- 51319  Rev.  1. 

9.  F.  R.  Gilmore  (1955),  Equilibrium  Composition  and  Thermo¬ 
dynamic  Properties  of  Air  to  2 4000^ K,  Rand  Report  RM-1543. 
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simply  as  a  flag  to  discriminate  which  cells  contained  air  and 
which  cells  contained  HE.  Pressures  in  mixed  cells  (w.iere 
0.1<f<0.9)  were  blended  according  to  Dalton's  law.  For  the 
present  study,  heat  release  due  to  turbulent  mixing  near  the 
HE/air  interface  was  neglected  (Sa=0) . 

c.  Results — As  a  baseline,  we  first  performed  an  inviscid 
calculation  of  this  HE-driven  blast  wave  (i.e.,  k=e=Ui^0). 

We  found  it  necessary  to  use  very  fine  zoning  for  this  problem. 
The  grid  spacing  was  initialized  as  one  tenth  of  one  percent  of 
the  shock  radius  (A/Rg=.001)  in  the  outer  region  of  the  charge 
(0. 08£r/Rsi.l>  0)  yielding  200  cells,  and  then  allowed  to 
increase  according  to  an  interest  rate  formula  (0<r/'R,q<0.8) 
over  50  cells,  with  a  grand  total  of  250  cells  inside  the 
charge.  The  cell  containing  the  primary  or  outer  shock  front 
was  moved  with  a  velocity  about  equal  to  the  instantaneous 
shock  velocity,  ws,  so  that  this  shock  remained  approximately 
fixed  in  the  grid.  Cell  interfaces  inside  the  shock  were 
assigned  velocities  linearly  proportional  to  the  radius: 

(Vi+i*  ”  ws  *  ri+>/RS  (11) 

Note  that  this  grid  law  produces  a  linear  dilation  of  the 
grid,  and  hence  preserves  the  initial  mesh  space  distribu¬ 
tion — a  very  useful  attribute  for  this  case.  It  keeps  the 
cells  packed  where  all  the  action  is — near  the  contact  surface 
and  shock  front! 

Figures  1-8  show  the  inviscid  flow  field  spatial  distri¬ 
butions  as  calculated  for  this  HE-driven  blast  wave  at  various 
instants  in  time,  starting  from  the  CJ  detonation  flow  field 
(cycle  0)  out  at  a  point  in  time  where  the  air  shock  has 
reached  about  two  charge  radii  (cycle  800).  Included  near 
the  bottom  of  each  figure  is  the  mesh  distribution  for  the 
last  displayed  time.  Note  in  particular  the  density  distri¬ 
butions  of  Figure  1  which  show  a  primary  shock  (Si) ,  followed 
by  a  contact  surface  (CS) ,  and  then  a  secondary  shock  (S2). 

The  latter  is  an  inward- facing  shock  that  is  being  swept 
outward  by  the  supersonic  flow.  At  later  times,  this 
secondary  shock  starts  to  move  inward  and  eventually  implodes. 

Next,  a  turbulence  calculation  was  performed.  Initial 
conditions  were  taken  as  the  flow  field  at  the  end  on  the 
aforementioned  inviscid  calculation  for  the  following 
reasons: 

•  At  this  point  the  three  discontinuities  had 
separated  from  each  other  and  were  well 
resolved  on  the  mesh. 

•  This  allowed  an  accurate  calculation  of  the 
pressure  and  density  gradients  which  fed  into 
the  turbulence  source  term,  G. 
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Figure  1.  Spatial  distribution  of  the  density  field  at  various  times 
(inviaeid  FCT  calculation  of  a  1-xb  spherical  PBX-9404 
charge  detonated  in  air) - 
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Figure  3.  Spatial  distribution  of  the  total  energy  field  at  various  times 
(inviscid  FCT  calculation  of  a  1-lb  spherical  PBX-9404 
charge  detonated  in  air) . 


The  values  of  the  k-e  model  parameters  were  taken  as  those 
established  from  turbulent  shear  flow 

Ci  =  1.44 

C2  =  1.92 

°T,8  =  1*0  for  $  =  h,k,f 

^T  > £  ~  1*3 

while  the  viscosity  coefficient,  Cy,  and  the  Rayleigh-Taylor 
coefficient,  C3,  were  assumed  to  be: 

Cy  =  0.009 

C3  =  1.0 

Figures  9  through  19  give  the  evolution  of  the  flow  field 
at  100  and  200  cycles  after  the  turbulence  terms  were  turned 
on.  These  figures  show  a  narrow  region  of  the  flow  field  in 
the  vicinity  of  the  contact  surface.  Consider  in  particular 
Figures  17,  18  and  19  which  depict  the  calculated  turbulent 
kinetic  energy,  k,  dissipation  rate,  e,  and  total  viscosity, 
Ut*  Comparing  these  parameters  with  the  density  distribution 
shows  that  the  turbulent  kinetic  energy  grows  rapidly  in 
the  vicinity  of  the  contact  surface  and  nowhere  else.  After 
200  cycles,  k^5  x  10®  erg/g  while  internal  and  kinetic 
energies  are  greater  than  I0l0  erg/g.  Note  that  the  dissi¬ 
pation  rate  is  also  growing?  hence,  the  evolution  of  the 
turbulence  kinetic  energy  will  depend  on  the  delicate  balance 
between  the  generation  terms  and  the  decay  terms.  These,  in 
turn,  depend  on  the  particular  values  chosen  for  C3  and  Cy. 

A  detailed  investigation  of  these  effects  is  an  appropriate 
subject  for  future  studies. 
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Pigure  13.  Spatial  distribution  of  the  pressure  field  at  various  times 
(continuation  cf  the  FCT  calculation  of  Figures  1-8  with 
the  k-e  turbulence  -nodel  activated)  . 


Figure  14.  Spatial  distribution  of  the  internal  energy  field  at 
time*  (continuation  of  the  FCT  calculation  of  Figure: 
with  the  k-e  turbulence  model  activated) . 


Spatial  distribution  of  the  species  fraction  field  at  various 
times  {continuation  of  the  FCT  calculation  of  Figures  1-8 
with  the  fe-E  turbulence  model  activated) . 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  shown  here  that  the  present  k-e  turbulence  model 
does  indeed  predict  the  rapid  growth  of  turbulence  in  the 
vicinity  of  the  contact  surface  for  HE-driven  blast  waves. 
Parametric  studies  should  be  performed  with  this  model  to 
pin  down  the  constants,  especially  C3  and  c^.  Such  kinematic 
calculations  should  be  verified  with  experimental  data.  In 
addition,  a  heat  release  model  should  be  formulated  and 
incorporated  into  the  hydrocode.  Calculations  including  such 
energetics  (i.e.,  afterburning)  should  then  be  performed  for 
explosives  which  afterburn  more  strongly  (e.g.,  TNT  and 
Pentolite)  and  more  weakly  (e.g.,  PBX-9404),  and  such  calcu¬ 
lations  should  then  be  compared  with  appropriate  experimental 
data  (pressure  records,  photography  and  thermal  radiation 
measurements  of  the  fireball)  and  inviscid  calculations  to 
determine  the  importance  of  afterburning  for  various  types 
of  explosives. 
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